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Abstract 


This  report  describes  some  of  the  problems,  achievements,  and  directions  of  investi¬ 
gation  of  three  areas  of  research.  The  first  is  the  investigation  of  parametric  nonlinear 
programming  problems  using  numerical  bifurcation  and  continuation  methods  and  with 
applications  to  design  optimization  and  parametric  control  systems.  The  second  part  cen¬ 
ters  on  investigation  of  various  numerical  methods  for  the  solution  of  nonlinear  optimal 
control  problems.  The  analysis  of  convergence  in  infinite  dimensional  spaces,  discretiza¬ 
tions,  and  numerical  implementations  are  in  progress  for  Newton’s,  penalty,  augmented 
Lagrangian,  and  interior  point  methods.  The  third  part  ot  this  research  program  is  the 
development  of  combinatorial  optimization  techniques  to  solve  the  central  problem  of  multi¬ 
target  tracking,  i.e.,  the  data  association  problem  of  partitioning  observations  into  tracks 
and  false  alarms.  The  problem  formulation,  algorithm  design,  and  real-time  solution  in¬ 
volve  techniques  from  probability  and  information  theory,  system  identification,  filtering, 
control  systems,  combinatorial  optimization,  and  advanced  computer  architectures,  includ¬ 
ing  massively  parallel  computers.  The  data  association  problem  for  general  multi-target 
tracking  problems  is  posed  as  a  class  of  multi-dimensional  assignment  problems.  The 
algorithms  under  development  are  based  on  a  recursive  Lagrangean  relaxation  scheme, 
construct  high  quality  suboptimal  solutions  in  real-time,  and  use  a  variety  of  techniques 
ranging  from  two  dimensional  assignment  algorithms,  a  conjugate  subgradient  method  for 
the  nonsmooth  optimization,  graph  theoretic  properties  for  problem  decomposition,  and  a 
branch  and  bound  technique  for  small  solution  components.  These  algorithms  are  being 
implemented  on  massively  parallel  computer  architectures  for  increased  performance  and 
real-time  identification. 
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1  Introduction 


This  report  describes  some  of  the  problems,  achievements,  and  directions  of  investi¬ 
gation  of  three  areas  of  research.  The  first  is  the  investigation  of  parametric  nonlinear 
programming  problems  using  numerical  bifurcation  and  continuation  methods  with  appli¬ 
cations  to  design  optimization  and  parametric  control  systems,  and  represents  a  potential 
for  a  real  extension  of  our  understanding  of  basic  phenomena,  global  sensitivity,  robustness, 
and  multiplicity  of  solutions  in  much  the  same  way  that  these  theoretical  and  numerical 
techniques  have  helped  the  understanding  of  dynamical  systems  and  nonlinear  equations. 
Thus  the  objective  in  this  aspect  of  the  research  program  is  to  develop  the  analytical  and 
numerical  techniques  to  map  out  regions  of  qualitatively  different  behavior  and  to  locate 
the  “stability”  boundaries  of  these  regions  in  parameter  space.  The  latter  is  important  be¬ 
cause  drastic  changes  in  the  optimum  occur  in  the  presence  of  singularities  which,  in  turn, 
define  these  “stability”  boundaries.  Such  knowledge  allows  for  the  uncertainty  in  system 
and  model  parameters  and  yields  information  about  the  expected  behavior  when  control 
parameters  are  varied  to  enhance  the  performance  of  the  system  under  consideration.  In 
addition  to  providing  a  global-like  sensitivity  analysis,  these  methods  are  quite  efficient  in 
computing  multiple  optima.  Several  model  problems  taken  from  the  very  active  area  of 
design  optimization  are  being  investigated  to  test  and  illustrate  the  value  and  applicability 
of  these  continuation  and  bifurcation  methods,  as  well  as  to  provide  motivation  and  focus 
for  further  development.  A  preliminary  theory  and  numerical  implementation  have  been 
completed  as  described  in  detail  in  Sections  2  and  3. 

The  second  part  centers  on  investigation  of  various  numerical  methods  for  the  solu¬ 
tion  of  nonlinear  optimal  control  problems.  The  analysis  of  convergence  in  infinite  dimen¬ 
sional  spaces,  discretizations,  and  numerical  implementations  are  in  progress  for  Newton’s, 
penalty,  augmented  Lagrangian,  and  interior  point  methods.  A  longer  term  goal  is  the  in¬ 
vestigation  of  parametric  problems  in  nonlinear  control  systems  including  but  not  limited 
to  the  nonlinear  optimal  control  problem.  Some  of  the  initial  results  in  this  direction  are 
described  in  Section  4. 

The  third  part  ot  this  research  program  is  the  development  of  combinatorial  op¬ 
timization  techniques  to  solve  the  central  problem  of  multi-target  tracking,  i.e.,  the  data 
association  problem  of  partitioning  observations  int  o  tracks  and  false  alarms.  The  problem 
formulation,  algorithm  design,  and  real-time  solution  involve  techniques  from  probability 
and  information  theory,  system  identification,  filtering,  control  systems,  combinatorial  op¬ 
timization,  and  advanced  computer  architectures,  including  massively  parallel  computers. 
The  data  association  problem  for  general  multi-target,  tracking  problems  is  formulated 
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as  a  class  of  multi-dimensional  assignment  problems.  The  algorithms  under  development 
are  based  on  a  recursive  Lagrangean  relaxation  scheme,  construct  high  quality  subopti- 
mal  solutions  in  real-time,  and  use  a  variety  of  techniques  ranging  from  two  dimensional 
assignment  algorithms,  a  conjugate  subgradient  method  for  the  nonsmooth  optimization, 
graph  theoretic  properties  for  problem  decomposition,  and  a  branch  and  bound  technique 
for  small  solution  components.  These  algorithms  are  being  implemented  on  massively  par¬ 
allel  computer  architectures  for  increased  performance  and  real-time  identification.  The 
current  status  and  results  of  this  research  effort  are  described  in  Section  5 

The  technical  information  concerning  publications,  lectures,  graduate  students,  col¬ 
leagues,  and  awards  is  contained  in  Section  7. 

2  Predictor-Corrector  Continuation  Algorithms 

Predictor-corrector  continuation  methods  for  tracing  solution  paths  of  an  under  de¬ 
termined  nonlinear  system  F(io)  =  0  where  F  :  lRn+1  — ►  IR"  have  proven  to  be  robust  and 
effective  procedures,  particularly  for  solving  problems  ranging  from  homotopy  methods  for 
nonlinear  equations,  continuum  mechanics,  and  optimization  to  the  study  of  parametric 
dependencies  in  dynamical  systems.  (The  parametric  optimization  problems  discussed  in 
the  next  section  represent  a  primary  source  of  problems  in  this  work.)  Although  these 
methods  are  quite  robust,  they  have  been  considered  slow  and  computationally  inten¬ 
sive,  primarily  because  of  the  extensive  linear  algebra  required  in  both  the  prediction  and 
the  correction  phases.  To  increase  the  computational  efficiency  and  robustness,  a  general 
Adams-Bashforth  predictor-corrector  continuation  procedure,  valid  for  both  homotopy  and 
parametric  problems,  has  been  developed  by  Lundberg  and  Poore  [6].  This  nonstandard 
ordinary  differential  equations  technique  employs  a  Newton-like  procedure  in  the  correc¬ 
tion  process  and  a  variable  order  and  an  adaptive  stepsize  control  in  the  prediction  phase 
to  efficiently  and  robustly  trace  out  paths  of  solutions  in  these  parametric  problems.  We 
[6]  have  sliowm  these  procedures  to  be  highly  efficient  and  in  many  respects  superior  to 
existing  path  following  algorithms. 

3  Parametric  Problems  in  Nonlinear  Programming  and  Control 

3A  Problem  Statement.  The  parametric  nonlinear  programming  problem  is  that  of 
determining  the  behavior  of  the  solution(s)  as  a  parameter  or  vector  of  parameters  cv  €  lRr 
varies  over  a  region  of  interest  for  the  problem 

Minimize 

Subject  To  h{  r ,  n  )  =  0  (3.1) 

r/(.r,o)  >  0 
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where  /  :  Itn+r  — ►  H,  Hu+r  — ►  1R9  and  g  :  Rn+r  — >  Rp  are  assumed  to  be  at  least  twice 
continuously  differentiable.  Most  applied  problems  contain  parameters;  some  may  be  fixed 
but  not  known  precisely  and  others  may  be  varied  to  enhance  the  performance  of  the 
system.  Good  local  information  about  the  rates  of  change  of  the  variables  with  respect  to 
the  parameters  can  be  obtained  by  differentiating  the  Karush-Kuhn-Tucker  conditions  and 
this  procedure  can  be  rigorously  justified  at  regular  points  of  (3.1)  by  the  implicit  function 
theorem.  Sensitivity  is  due  to  a  “nearby”  singularity  in  the  problem,  and  thus  we  have 
investigated  these  singularities  extensively  [10,18,21,22].  These  singularities  arise  from  a 
loss  of  strict  complementarity,  a  loss  of  the  linear  independence  constraint  qualification,  or 
the  singularity  of  Hessian  of  the  Lagrangian  on  the  tangent  space  to  the  active  constraints 
[18], 

To  our  knowledge,  numerical  continuation  and  bifurcation  techniques  have  not  been 
systematically  developed  for  the  fully  constrained  problem.  These  same  methods,  however, 
have  been  highly  successful  in  the  numerical  study  of  dynamical  systems  and  nonlinear' 
equations,  and  thus  the  current  objective  is  to  develop  numerical  continuation  and  bifur¬ 
cation  techniques  in  parametric  nonlinear  programming  in  the  near  term  and  in  control 
systems  in  the  longer  term.  In  addition  to  yielding  a  global-like  sensitivity  analysis  of  the 
parametric  problem,  numerical  continuation  and  bifurcation  methods  also  yield  an  effective 
method  for  computing  multiple  solutions. 

To  test  and  illustrate  the  value  and  applicability  of  our  methods  as  well  as  to  provide 
motivation  and  focus  for  further  development,  we  are  investigating  several  model  prob¬ 
lems  taken  from  the  very  active  area  of  design  optimization.  We  are  confident  that  the 
techniques  currently  under  development  for  parametric  constrained  optimization  problems 
will  be  very  useful  for  investigating  sensitivity,  stability,  and  multiple  optima  in  structural 
design  and  control  problems.  Our  first  publication  on  this  problem  is  that  of  Lundberg 
and  Poore  [7]. 

3B  Status  of  the  Numerical  Algorithms.  To  develop  a  numerical  continuation 
and  bifurcation  approact  to  the  parametric  nonlinear  programming  problem,  a  theoretical 
development  of  the  singularities  in  parametric  nonlinear  programming  and  a  thorough 
understanding  of  a  good  continuation  code  are  required.  Tiahrt  and  Poore  [10,18,21,22] 
have  investigated  singularities  and  persistence  of  the  minima  as  well  as  critical  point  type  in 
nonlinear  parametric  programming  and  Lundberg  and  Poore  [6]  have  developed  an  Adams- 
Bashfort.h  predictor-corrector  continuation  method  which  employs  variable  order  and  an 
adaptive  stepsize  control  to  efficiently  and  robustly  trace  out  paths  of  solutions  in  these 
parametric  problems.  (This  nonstandard  ODE  technique  uses  a  Newton-like  method  in 
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the  correction  process.)  To  explain  the  current  status  of  these  numerical  investigations,  we 
first  pose  the  parametric  programming  problem  as  a  system  of  nonlinear  equation^  [18,21]. 
From  the  Fritz  John  first-order  necessary  conditions,  there  exist  q  +  p  +  1  real  numbers  in 
v,  A  =  (Aj, . . . ,  A q)  and  p  =  (p\, . . .  ,pp),  not  all  zero,  such  that 


F(x,v,\,p]a) 


'  Vx£(:r,i/,  \,p\a) 
—h(x,  a) 
Mg(x,  a) 

.v2  +  pTp  -|-  ArA  -  /?o 


(3.2) 


where  £  =  vf(x,<y)—  <  A ,h(x,a)  >  —  <  p,g(x,a)  >  is  the  Lagrangian.  In  the  presence 
of  a  constraint  qualification,  the  usual  normalization  is  to  delete  the  last  equation  and 
set  v  =  1;  however,  the  loss  of  the  linear  independence  constraint  qualification  leads  to 
a  singularity  which  is  generally,  but  not  always,  characterized  by  multipliers  tending  to 
infinity.  This  nonstandard  normalization  for  the  multipliers  n,  A,  and  p  replaces  multipliers 
tending  to  infinity  by  u  — ►  0.  This  system  system  contains  all  solutions  of  the  Fritz  John  or 
Karush-Kuhn-Tucker  first  order  necessary  conditions  as  well  as  minima,  maxima,  saddle 
points,  and  both  feasible  and  infeasible  solutions. 

To  explain  the  development,  please  keep  in  mind  that  singularities  in  the  above  system 
can  only  arise  from  a  loss  of  strict  complementarity,  a  loss  of  the  linear  independence 
constraint  qualification,  or  the  singularity  of  Hessian  of  the  Lagrangian  on  the  tangent 
space  to  the  active  constraints  [18,21].  The  easiest  singularity  to  treat  is  that  of  the  loss 
of  strict  complementarity,  which  corresponds  to  a  bifurcation  point  whose  branches  can  be 
delineated  by  activating  and  deactivating  the  particular  constraint(s)  in  question.  Since 
the  signs  of  the  inequality  constraints  and  corresponding  multipliers  can  be  monitored  to 
detect  these  bifurcations,  the  problem  of  path  following  and  bifurcation  detection  in  one 
parameter  can  be  simplified  to 


F(x,i',\,p;a)  = 


Vj£(r,  u.  A,  p\ «) 

—h(x,n) 

-g(x,a) 

.  u‘2+  <  //,  p  >  +  <  A,  A  >  -fil . 


(3.3) 


where  g  and  p  are  obtained  from  g  and  p  by  deleting  constraints  and  multipliers  cor¬ 
responding  to  inactive  constraints.  (This  is  essentially  an  “active  set  strategy.")  Given  a 
general  predictor-corrector  continuation  method  [6],  the  development  of  numerical  continu¬ 
ation  and  bifurcation  techniques  for  constrained  optimization  thus  begins  with  the  tailoring 
of  the  numerical  linear  algebra  and  singularity  detection  methods  to  this  formulation  of 
the  problem. 
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In  the  continuation  process,  one  also  adds  an  augmenting  equation  to  specify  how  the 
correction  is  made  back  to  the  path  after  a  prediction,  so  that  the  coefficient  matrix  in  the 
linear  systems  that  must  be  solved  at  each  point  along  the  path  are  of  the  form 

,  \W  B]  .  ...  [V*£  -A]  .... 

J  =  cT  D  where  W=  -At  0  (3-4) 

and  At  represents  the  derivative  with  respect  to  x  of  the  active  constraints.  The  last  two 
rows  and  columns  in  J,  represented  by  the  matrices  CT ,  D,  and  B,  arise  from  the  multiplier 
normalization  above  and  the  augmenting  equation.  An  effective  framework  for  developing 
both  the  numerical  linear  algebra  and  singularity  detection  is  a  slight  modification  of 
Keller’s  bordering  algorithm  [4,5],  which  reduces  the  linear  systems  involving  J  to  those 
involving  W.  Then  null  and  range  space  methods  and  the  direct  factorization  of  W,  which 
are  extensively  developed  and  used  in  optimization  [1,  Section  10.2],  can  be  adapted  to  the 
continuation  problem  with  one  modification.  Since  in  the  continuation  process  one  follows 
paths  of  all  types  of  critical  points  including  minima,  the  Hessian  of  the  Lagrangian  on 
the  tangent  space  to  the  active  constraints  need  no  longer  be  positive  definite.  Thus  the 
linear  algebra  methods  need  to  be  modified,  e.g.  by  replacing  a  Cholesky  factorization  by 
a  LDLT  factorization  using  the  Bunch-Kaufman  algorithm  [2]. 

As  discussed  above  the  loss  of  strict  complementarity  can  be  detected  by  monitoring- 
sign  changes  in  an  inactive  inequality  constraints  and  multipliers  corresponding  to  active 
equality  constraints.  Thus  methods  must  be  developed  for  the  remaining  two  singularities. 
Changes  in  the  signature  of  the  Hessian  of  the  Lagrangian  on  the  tangent  space  to  the 
active  constraints  can  be  accomplished  directly  for  null  space  methods,  and  using  results 
on  Schur  complements  [19]  and  closely  related  inertia  results  [3],  we  [8]  have  been  able  to 
develop  efficient  methods  for  this  detection  for  both  the  range  space  method  and  the  direct 
factorization  of  W.  Loss  of  the  linear  independence  constraint  qualification  occurs  when 
v  — *  0  since  then  there  is  a  nonzero  vector  A  such  that  AX  =  0.  We  are  also  investigating 
methods  based  on  the  factorization  of  the  constraint  matrix  AT .  The  systematic  devel¬ 
opment  of  these  methods  will  be  reported  in  a  forthcoming  work  of  Lundberg  and  Poore 

[7,8]. 

3C  A  Numerical  Example  from  Design  Optimization.  As  a  simple  illustration  of 
the  above  procedures  we  consider  the  following  problem  from  design  optimization  [20]: 

Minimize  <7 

Subject  To  V,iE(d,  h;p)  =  0  (3.5) 

0  <  h  <  1.5 


i 
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where  E(d,h;p)  =  -pd  +  (Vl  +  h2  -  +  (h  -  d )2j  /\/l  +  h2 ,  and  p  is  the  parameter. 

This  problem  is  used  to  model  the  determination  of  the  unloaded  height  h  of  a  simple 
two  bar  planar  truss  with  semi-span  1,  which  minimizes  the  displacement  d  under  a  fixed 
load  p.  The  solution  paths  z(p)  =  (d(p),  h(p),  A(p),  p(p))  of  (3.5)  were  tracked  using  our 
continuation  method.  The  following  plot  gives  the  displacement  d  as  p  varies. 


Fig  3.1:  Feasible  Solutions  of  (3.5) 
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This  plot  represents  a  projection  of  the  feasible  solutions  of  (3.5)  into  the  (p,  d)  plane, 
and  the  dot  labeled  with  e,  /  and  g  indicates  three  distinct  bifurcation  points  with  p  =  d  =  0 
and  h  =  0,  1.41  and  1.5,  respectively.  Bifurcation  points  a,c,  e  and  g  result  from  a  loss 
of  strict  complementarity  in  which  an  inequality  constraint  becomes  weakly  active.  The 
path  of  maximizers  branching  from  point  a  corresponds  to  h  <  1.5  active  and  pi  <  0.  and 
changes  type  at  the  singular  point  g,  becoming  the  path  of  minimizers  labeled  G.  The 
other  path  branching  from  point  a  passes  through  /,  across  which  the  one  eigenvalue  of 
changes  sign.  At  /  there  is  a  change  in  type  resulting  in  the  path  of  maximizers 
labeled  F. 

Extreme  sensitivity  of  the  solution  of  (3.5)  to  variations  in  p  occurs  at  the  fold  points 
b  and  d  (p  =  ±.3704)  at  which  there  is  a  loss  of  linear  independence  in  the  active  constraint 
gradients,  and  p2  = *0-  This  is  also  the  case  at  the  bifurcation  point  e,  where  in  addition, 
strict  complementarily  is  violated.  One  cannot  compute  near  or  past  these  points  without 
the  normalization  v2  +  \T \  +  pT p  —  do  =  0,  since  near  these  points  an  unnormalized 
multiplier  is  unbounded.  When  the  system  is  at  a  state  near  these  points,  small  variations 
in  load  p  can  result  in  very  large  changes  in  the  solution,  or  the  loss  of  (local)  existence 
of  a  solution.  The  latter  case  is  illustrated  near  b  where  increasing  the  parameter  past 
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p  =  .3704  results  in  the  loss  of  the  solution  and  a  “snap  through1'  of  the  truss  to  a  state 
represented  by  path  D.  Similar  behavior  occurs  near  d  and  e. 

Not  pictured  above  are  branches  of  infeasible  solutions  of  (3.5)  emerging  at  a,c,  e 
and  g  ( h  >  1.5  or  h  <  0).  (In  some  problems  such  paths  may  provide  the  opportunity 
for  further  branching  to  other  feasible  paths.)  A  path  of  feasible  singular  points  with 
p  =  d  =  0  branches  from  e  through  /  to  g ,  and  can  be  parameterized  by  p2-  The  solutions 
to  (4.1)  need  only  be  stationary  points  of  the  potential  energy  E(d,h;p).  However,  all 
path  segments,  exclusive  of  the  segments  from  d  to  e  and  from  c  to  b,  do  correspond  to 
physical  states  of  the  system  (where  E  is  minimized). 

Finally,  note  the  multiple  solutions  in  the  diagram  which  are  easily  computed  via  these 
continuation  procedures. 

4  Nonlinear  Optimal  Control 

The  optimal  control  problem  under  consideration  in  this  work  can  be  described  as 
Minimize  J[x,tt]  :=  tp(x(t0),x(tl))  +  f  f0(t,x,u)dt 

Ji  0 

Subject  To  x  —  f(t,x(t),a(t)) 

B(x(to),.r(tl))  =  0 
h(t,x,  u)  =  0 
<7(t,.r,  u)  >  0 
u  e  Q 

(x,u)  e  trl-°°([*o,*i],Rn)  x  T°°([*o,*I],Rm) 

where  x  is  an  ’-vector,  u  is  an  m-vector,  B  is  a  boundary  operator,  Q  is  a  closed  convex 
set,  and  VF1,p([t0,  t\  ],  1R"  )  is  the  usual  Sobolev  space  which  can  be  characterized  via  the 
Sobolev  imbedding  theorem  as  consisting  of  those  absolutely  continuous  vector  functions 
with  the  first  derivative  in  T;'([/y,  / 1 ] ,  IR").  The  functions  ^?,/o,/,  B,  h,  and  g  are  assumed 
to  be  at  least  C2  with  respect  to  their  arguments.  In  this  formulation,  we  assume  that 
the  t o  and  t.  \  are  fixed;  however,  we  stress  that  the  more  general  problem  in  which  <p  = 
■p(to,  x{t0 ),  b ,  x(t j ))  and  the  end  points  (b,  x(t, ))  are  allowed  to  vary  can  be  transformed 
into  this  problem  by  a  standard  transformation  that  introduces  two  new  state  variab.  ’. 

Our  interest  in  this  problem  is  two  fold.  First,  working  with  W.  W.  Hager  of  the 
University  of  Florida,  we  are  investigating  the  convergence  of  various  numerical  methods 
(Newton’s,  penalty,  augmented  Lagrangian,  interior  point  methods)  in  the  appropriate 
infinite  dimensional  spaces  and  will  then  work  on  the  discretization  and  numerical  solution 
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of  these  problems.  One  paper  has  been  submitted  for  publication  [16],  one  is  in  preparation 
[17],  and  a  third  is  planned 

The  second  area  of  interest  is  in  the  parametric  problem  obtained  by  inserting  pa¬ 
rameters  in  the  above  problem.  The  interaction  of  multiple  and  bifurcating  states  in  the 
absence  of  controls,  periodic  phenomena,  chaotic  behavior,  and  bifurcating  controls  arising 
from  the  dynamical  systems  and  holonomic  constraints  is  open  to  investigation.  Given  a 
certain  phenomena  arising  from  a  dynamical  system,  the  problem  may  be  to  control  this 
phenomena,  to  determine  multiple  solutions,  or  to  investigate  the  dependence  of  a  solution 
on  the  system  parameters  over  a  wide  range,  i.e.  global  sensitivity.  (The  latter  is  also 
important  in  adaptive  control.)  The  development  and  use  of  theoretical  and  numerical 
bifurcation  and  continuation  methods  in  dynamical  systems  and  nonlinear  equations  has 
been  spectacularly  successful  in  analyzing  and  understanding  the  phenomena  represented 
by  these  systems,  but  we  know  ol  no  systematic  treatment,  or  works  on  the  constrained 
nonlinear  parametric  control  problem  parallelling  that  found  in  dynamical  systems.  Thus 
a  long  term  goal  of  this  research  program  will  be  the  investigation  of  parametric  problems 
in  nonlinear  control  systems  including  but  not  limited  to  the  nonlinear  optimal  control 
problem. 

5  Combinatorial  Optimization  and  Multi-Target  Tracking 
5A  Problem  Statement. 

The  third  part  ot  this  research  program  is  the  development  of  combinatorial  opti¬ 
mization  techniques  to  solve  the  central  problem  of  multi-target  tracking,  i.e.,  the  data 
association  problem  of  partitioning  observations  into  tracks  and  false  alarms.  Although 
combinatorial  optimization  is  the  natural  framework  for  the  formulation  of  these  problems, 
the  corresponding  techniques  have  long  been  considered  computationally  too  intensive  for 
real-time  applications  and  for  good  reason.  The  resulting  optimization  problems,  which 
are  formulated  here  as  inulti-dimensional  assignment  problems,  are  NP-hard.  To  further 
appreciate  the  difficulties,  one  only  has  to  examine  the  trade-offs  between  two  current 
methods  in  multiple  target  tracking:  track  while  scan  and  batch.  For  the  former,  one 
essentially  extends  tracks  a  scan  at  a  time  using  for  example  a  two  dimensional  assignment 
or  a  greedy  algorithm.  This  methodology  is  real-time,  but  results  in  a  large  number  of 
partial  and  incorrect  assignments,  and  thus  incorrect  track  identification.  The  fundamen¬ 
tal  difficulty  with  this  approach  is  flier e  is  simply  not  enough  information  in  one  scan  at 
a  time  processing  to  properly  partition  the  observations  into  tracks  and  false  alarms.  To 
obtain  the  required  information,  one  needs  to  consider  several  scans  all  at  once,  i.e.  the 
batch  approach,  but  it  is  this  batch  approach  that  is  considered  computationally  too  int.cn- 
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sive  for  practical  real-time  applications.  Given  the  ever-present  need  to  identify  all  tracks 
in  real-time ,  the  challenge  to  combinatorial  optimization  is  to  design  fast  algorithms  for 
advanced  computer  architectures  that  will  solve  the  underlying  data  association  problem, 
and  thus  the  track  identification,  in  real-time. 

What  we  have  accomplished  is  the  real-time  or  near  real-time  of  many  data  association 
problems  by  exploiting  sparsity  and  problem  decomposition,  by  developing  a  Lagrangian 
relaxation  algorithm  for  the  construction  of  ‘high  quality1  suboptimal  solutions,  and  by 
using  advanced  computer  architectures  such  as  massively  parallel  architectures.  In  the 
subsections  to  follow,  the  problem  formulation,  algorithms,  current  results,  and  parallel 
computing  issues  are  discussed. 

5B  Physical  Model.  A  number  of  objects  or  targets,  say  100  -  10,000,  are  assumed  to 
obey  some  underlying  physics  from  which  one  can  formulate  physical  Lavs  for  the  equations 
of  motion.  Most  models  for  real-time  application  arc  based  on  linear  systems  theory 
wherein  each  target  is  assumed  to  obey  an  equation  of  the  form 

=  F(t)x  +  T(t)u  +  G(t)w 


z  —  H(t).v  +  v 

where  F(t).  T( / ),  G(t),  H(t)  are  assumed  to  be  known,  x  is  the  state  variable,  u  is  the 
input  or  control  function,  iv  is  the  input  or  process  noise,  v  is  the  observation  error,  and  c 
is  the  output  or  observed  quantity.  What  differentiates  one  target  from  another  might  be 
initial  state  values  or  parameters  in  F{t).  Tft),  G(t),  H(t );  these  time  dependent  matrices 
can  also  vary  from  target  to  target  in  a  more  general  way. 

At  a  set  of  scan  times  {  f  t ,  pictures  of  the  objects  are  taken,  and  the  observations 
are  recorded  as  {;*  } i  k=i)  for  scan  time  /*..  Here,  K  represents  the  number  of  scans.)  Ar*  is 
the  number  of  observations  on  scan  k ,  and  the  zero  index,  7*  =  0,  corresponds  to  a  dummy 
or  missed  observation. 

Given  this  description,  the  data  association  problem  is  that  of  partitioning  the  obser¬ 
vations  into  tracks  and  false  alarms  [11,12]  in  such  a  way  that  the  paths  or  tracks  of  the 
objects  can  be  identified.  Smoothing,  filtering,  and  prediction  techniques  are  then  used  to 
obtain  further  information  about  past.,  present,  and  future  states  of  the  objects. 

5C  Mathematical  Formulation  of  the  Data  Association  Problem.  In  what,  fol¬ 
lows,  the  term  track  of  observations  is  used  to  denote  a  sequence  of  observations  { 
one  from  each  scan,  that  might  be  generated  by  the  target  or  object.  Since  the  potential 
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number  of  tracks  of  observations  is  too  large  to  consider  computationally,  a  gating  pro¬ 
cedure  is  employed  to  remove  unlikely  tracks  of  observations  and  thus  introduce  sparsity 
into  the  problem  [11]. 

Given  a  track  of  observations  { }  £'_  j ,  the  filtering  problem  is  to  estimate  the  state 
up  to  the  current  time  f  *.  The  first  use  of  filtering  and  system  identification  is  in  the  devel¬ 
opment  of  a  score  function  since  measurement  error  in  the  track  of  observations  is  scored 
against  the  filtered  track.  Specifically,  Kalman  filtering,  recursive  system  identification, 
and  adaptive  filtering  techniques  are  particularly  relevant  to  this  problem. 

Using  K  scans  of  information,  the  next  task  is  to  formulate  a  I\  -dimensional  assign¬ 
ment  problem  whose  solution  gives  an  optimal  partitioning  of  observations  into  tracks  and 
false  alarms.  Given  a  track  of  observations  (z^ , . . . ,  z^r ),  define  the  0-1  variable 

_  [  1  if  , . . . ,  z\yK  )  is  assigned  to  a  track, 
s  [^  0  otherwise. 

The  score  of  the  assignment  of  the  observations  (zj  , . . . ,  z-kk  )  to  a  track  is  defined  to  be 

. if  (**,,••■,*£)  Passes  gating, 

"''2 . 'K  l°o  if  (*',,-••,*£)  fails  gating, 

where  p,,  . ,K  is  a  composite  probability  density  function  involving  probabilities  for 

measurement  error,  false  alarms,  missed  detections,  track  life,  initiation,  and  finite  sensor 
resolution  [13]. 

The  problem  of  assigning  0  or  1  to  all  the  variables  Sj, . jK.  in  such  a  way  that  each 

actual  observation  is  assigned  to  exactly  one  track  total  cost  is  minimized  is  called  the 
assignment  problem  which  can  be  formulated  as 


A',  NK 

Minimize  r(  c )  =  ^  c„ ,  . ,K 

1 1=0  ijc=0 

A’j  Nk 

Subj.  To  X!  2^  r" . **■  =  =  1 . Nl' 
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An  important  feature  of  this  formulation  is  that  the  number  of  targets  is  not  required 
a  priori  and  is  determined  as  part  of  the  solution  process. 

5D  Algorithms  for  the  Construction  of  Real-Time  Solutions.  The  basic  scheme 
currently  used  employs  preprocessing  in  the  form  of  gating  and  problem  decomposition. 
Then  the  sparse  decomposed  problems  are  solved  by  a  recursive  Lagrangean  relaxation 
scheme.  A  A'-dimensional  assignment  problem  is  relaxed  to  a  (K  —  l)-dimensional  one 
by  incorporating  one  set  of  constraints  into  the  objective  function  using  a  Lagrangean 
relaxation  of  this  set.  Given  a  solution  of  the  (A  —  l)-dimensional  problem,  a  feasible 
solution  of  the  A -dimensional  problem  is  then  reconstructed.  The  ( K  —  l)-dimcnsional 
problem  is  solved  in  a  similar  manner  and  the  process  is  repeated  until  one  reaches  the  two- 
dimensional  problem  which  is  solved  exactly.  The  duality  gap  in  this  process  is  generally 
quite  small  and  one  obtains  in  general  an  e-optimal  solution.  The  full  technical  description 
can  be  found  in  the  forth  coming  paper  of  Poore  and  Rijavec  [11], 

5E  A  Case  Study.  The  following  three  tables  illustrate  the  solution  quality,  current 
timings,  and  expected  timings  for  the  identification  of  100  tracks  consisting  of  straight 
lines  in  two  dimensional  space-time.  The  two  tables  below  present  a  numerical  comparison 
for  the  straight  lines  whose  intercepts  and  slopes  are  uniformly  distributed  [0,1000]  and  [- 
0.2,0. 2],  respectively.  The  observation  errors  are  assumed  to  be  Gaussian  random  variables 
with  a  zero  mean  and  standard  deviation  a.  In  the  following  tables  the  maximum  error 
of  3 cr  is  defined  as  the  mean  distance  between  the  tracks  at  the  initial  time  and  ranges 
between  1%  and  50%.  The  scan  times  are  spaced  40  seconds  apart.  These  problems  are 
scale  invariant  for  rnmazAt  =  0.2  x  40  where  At  denotes  the  time  between  scans.  Thus 
if  the  observations  are  taken  every  5  seconds,  the  slopes  can  range  between  -1.6  and  1.6. 
All  computations  were  performed  on  Silicon  Graphics  Personal  IRIS,  and  twenty  (20)  test 
problems  were  ratrdomly  generated  and  solved  to  obtain  the  averaged  results  given  in  the 
two  tables  below.  Tints  an  identification  reading  of  99.95%  implies  that  all  but  one  track 
orrt  of  2000  was  correctly  identified. 


%  max.  error 

3  scans 

4  scans 

5  scans 

6  scans 

7  scans 

8  scans 

1.0 

100.00 

100.00 

100.00 

99.90 

100.00 

100.00 

5.0 

99.30 

99.60 

99.80 

99.85 

100.00 

99.95 

10.0 

97.10 

99.15 

99.80 

99.65 

99.85 

99.70 

20.0 

95.75 

98.20 

98.85 

98.85 

99.05 

99.25 

30.0 

94.50 

96.35 

97.90 

97.70 

98.80 

98.75 

40.0 

92.10 

94.20 

96.30 

97.50 

97.30 

99.25 

50.0 

93.50 

94.20 

95.60 

95.10 

97.25 

- 

Solution  Quality:  %  of  Tracks  Correctly  Identified 
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The  maximum  error  of  1%  corresponds  to  a  high  signal  to  noise  ratio,  whereas  50%  cor¬ 
responds  to  a  very  low  signal  to  noise  ratio.  The  identification  is  of  exceptional  quality 
over  all  signal  to  noise  ratios  and  gradually  improves  as  one  moves  across  the  table.  The 
reason  for  incomplete  identification  is  that  one  encounters  regions  of  high  contention  where 
many  tracks  cross.  This  difficulty  is  resolved  locally  on  subsequent  scans,  but  other  regions 
appear.  In  the  next  table  the  solution  times  are  presented. 


%  max.  error 

3  scans 

4  scans 

5  scans 

6  scans 

7  scans 

8  scans 

1.0 

0.05 

0.05 

0.06 

0.06 

0.07 

0.07 

5.0 

0.19 

0.31 

0.42 

0.59 

1.01 

1.61 

10.0 

0.66 

1.31 

1.94 

3.24 

5.01 

6.55 

20.0 

1.25 

3.32 

7.15 

16.60 

27.86 

53.45 

30.0 

1.96 

6.01 

18.10 

52.58 

118.15 

276.58 

40.0 

2.14 

8.50 

28.93 

90.69 

281.15 

506.11 

50.0 

2.37 

11.03 

42.98 

156.76 

379.62 

- 

Current  Solution  Times  in  Seconds 


We  have  mentioned  real-time  in  this  work;  let  us  now  be  specific.  Suppose  a  radar  sweep 
takes  5  to  10  seconds.  The  objective  then  is  to  process  as  many  scans  as  possible  between 
sweeps  to  improve  identification  and  solve  the  problem  in  the  alloted  time.  The  above 
table  gives  some  idea  of  the  capability  at  the  present  for  100  targets.  To  appreciate  these 
timings  further,  one  must  consider  the  several  factors  that  can  significantly  improve  the 
speed.  The  research  code  based  on  the  current  algorithms  is  designed  for  adaptability,  not 
speed.  (A  special  purpose  three  arid  four  dimensional  assignment  code  ran  between  six  and 
ten  times  faster  than  our  recursive  Iv-dimensional  code  used  to  generate  these  numbers.) 
Perhaps  the  most  significant  speedup  can  be  achieved  by  what  has  not  been  used.  Given 
a  specific  K  dimensional  assignment  problem,  one  most  often  has  a  suboptimal  or  optimal 
solution  of  a  closely  related  I\  or  K  —  1  dimensional  problem  and  such  a  solution  should 
be  used  as  a  hot  start.  For  example,  in  going  from  k  scans  of  observations  to  k  +  1  scans, 
almost  all  tracks  have  been  identified,  but  we  are  not  making  use  of  these  hot  starts.  The 
following  table  gives  the  expected  solution  times  for  this  problem. 
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%  max.  error 

3  scans 

4  scans 

5  scans 

6  scans 

7  scans 

8  scans 

1.0 

0.01 

0.01 

0.01 

0.01 

0.01 

0.01 

5.0 

0.02 

0.03 

0.04 

0.04 

0.05 

0.05 

10.0 

0.05 

0.07 

0.09 

0.09 

0.10 

0.15 

20.0 

0.07 

0.10 

0.15 

0.20 

0.20 

0.30 

30.0 

0.10 

0.15 

0.30 

0.50 

0.80 

1.20 

40.0 

0.10 

0.20 

0.40 

0.60 

1.40 

3.00 

50.0 

0.10 

0.20 

0.50 

1.20 

3.00 

4.00 

Expected  Solution  Times  in  Seconds 

5F  Parallel  Implementations.  The  multi-target  tracking  algorithms  are  in  the  pro¬ 
cess  of  being  parallelized  on  a  massively  parallel  SIMD  architecture.  The  goal  has  been  to 
investigate  and  understand  the  issues  and  algorithms  in  a  parallel  computing  environment. 
The  machine  chosen  for  algorithm  development  and  testing  was  the  Connection  Machine 
CM2  at  the  National  Center  for  Atmospheric  Research  (NCAR)  in  Boulder,  Colorado.  The 
Connection  Machine  consists  of  a  large  number  (up  to  64K)  of  1-bit  serial  processors,  each 
with  64K  bits  of  local  memory.  For  the  purpose  of  this  project,  only  a  small  number  of  the 
Connection  Machine  processors  were  used,  configured  into  a  linear  array  as  requested  by 
IBM  of  Ovvego,  New  York.  The  association  of  regions  of  space  to  processors  and  then  the 
load  balancing  for  efficiency  will  not  be  discussed,  since  problem  formulation  is  still  under 
development. 

Two  of  the  main  issues  in  implementing  algorithms  on  a  massively  parallel  issues  arc 
balancing  the  computational  loads  across  the  processors  and  minimizing  the  communica¬ 
tion  paths.  To  balance  computational  loads  we  have  developed  a  dynamic  load  balancing 
strategy  that  is  purely  local  in  nature,  i.e.  does  not  assume  the  existence  of  the  front 
end  computer.  Various  improvement  of  this  algorithm,  depending  on  a  particular  machine 
characteristics,  are  being  investigated.  Next,  a  comprehensive  strategy  for  fitting  the  track¬ 
ing  problem  to  a  massively  parallel  SIMD  architecture  has  been  developed.  This  strategy 
starts  with  dividing  the  target  space  into  slices  that  are  then  allocated  to  individual  pro¬ 
cessors.  This  strategy  is  consistent  with  the  objective  of  minimizing  communication.  As 
an  added  benefit,  the  computing  loads  arising  from  such  allocation  are  more  likely  to  be 
balanced,  making  the  load  balancing  algorithms  cheaper  to  execute. 

The  issue  of  parallel  indexing  (i.e.  indexing  parallel  arrays  by  parallel  indices)  has 
turned  out  to  be  fundamentally  important  for  efficient  implementation  of  tracking  algo¬ 
rithms  on  a  massively  parallel  SIMD  architecture.  Most  algorithms  require  at  least  some 
parallel  indexing,  while  some  are  very  heavily  dependent  upon  it.  The  architecture  that 
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does  not  provide  an  efficient  parallel  indexing  will  not  be  well  suited  for  these  tracking 
algorithms. 

Specifically  some  of  the  modules  that  have  been  implemented  todate  include: 

a  A  tracking  model  to  generate  the  observations.  This  model  is  implemented  in  parallel 
and  incorporates  targets  moving  with  constant  speeds  in  one  dimensional  space,  target 
initiation  and  termination  at  any  time,  probability  of  detection  of  less  then  one,  false 
alarms  (clutter)  and  finite  sensor  resolution, 
b  Problem  formulation  module  that  formulates  a  multi- dimensional  assignment  problem 
on  the  basis  of  observations  and  known  model  parameters. 
c  A  decomposition  algorithm  to  identify  disjoint  components  of  the  assignment  problem. 
d  A  branch  and  bound  algorithm  for  small  multi- dimensional  assignment  problems. 
e  The  auction  algorithm  of  Bertsekas  for  the  two  dimensional  assignment  problems. 
f  A  dynamic  load  balancing  algorithm. 
g  Various  sorting  and  utility  modules. 

Our  current  thinking  leads  to  the  conclusion  that  the  preferred  architecture  for  im¬ 
plementing  multi-target  tracking  algorithms  is  a  shared  memory  MIMD  architecture  with 
large  number  of  processors,  e.g.,  BBN  Butterfly.  However,  the  algorithms  are  also  well 
suited  to  implementation  on  a  massively  parallel  linear  SIMD  machine.  Such  machines  are 
much  less  complex  and  thus  likely  to  be  both  smaller  in  size  and  more  economical  to  build. 
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